function [nIters, errRel] = testBlockMG_node2(n, TOL, maxIt)
% A driver file to illustrate the usage of block-MG for Michael
%
% INPUT
%  n     : the number of cells, must be of the form 2^m.
%  TOL   : relative error threshold to stop the iteration.
%  maxIt : maximum number of iterations
%
% OUTPUT
%  nIters : the number of iterations performed
%  errRel : the relative error of the solution.
%
% example:
% [nIters, errRel] = testBlockMG_node(1024) should yield the following:
%   iter #1: max-norm of Residual=0.034663; max-norm of error=2.8076e-06
%   iter #2: max-norm of Residual=0.0011727; max-norm of error=4.3362e-08
%   iter #3: max-norm of Residual=4.1648e-05; max-norm of error=1.5591e-09
%   iter #4: max-norm of Residual=1.7096e-06; max-norm of error=6.2258e-11
%   iter #5: max-norm of Residual=7.0835e-08; max-norm of error=2.4709e-12


% This subroutine assembles the Laplacian operator with homogeneous Dirichlet BC.
    function Lap = LaplacianOperator(n, Dim, h)
        % return the 2nd-order Laplacian operator
        %   with homogenuous Dirichlet BC.
        % assuming cell-centered variables and the vector is columnwise
        
        L   = sparse([],[],[],n,n);
        nm1 = ones(n-1,1);
        L   = L + diag(-2*ones(n,1), 0) + diag(nm1, 1) + diag(nm1, -1);
        L = L/h/h;
        
        nP = n^Dim;
        if Dim==1
            Lap = L;
        elseif Dim==2
            In  = sparse(1:n, 1:n, ones(n,1), n, n);
            Lap = kron(L, In) + kron(In, L);
        else
            error('Unsupported value of Dim!');
        end
    end

if nargin<3
    maxIt = 30;
end
if nargin<2
    TOL = 1.e-6;
end

% the label for node-centering.
%  nVoffset==+1 : the boundary nodes are variables;
%  nVoffset==-1 : the boundary nodes are NOT variables.
nVoffset = +1;

addpath('./core');
global mg_matrix;
global mg_param;
mg_param.nu1 = 1;            % number of pre smooths
mg_param.nu2 = 1;            % number of post smooths
mg_param.minlevel = 1;       % level of coarsest grid
mg_param.bx  = 5;            % block size x           
mg_param.by  = 5;            % block size y
mg_param.sx  = 4;            % skip size x
mg_param.sy  = 4;            % skip size y
cycleType = 'f';             % MG cycle type.
                             % the other choice is 'v'
% rectangular grids with node-centered unknowns.
% n is the number of cells.
bds = [-2, 2];
h    = 1/n;
x1d  = linspace(bds(1), bds(2), n+1);
% the following takes the internal nodes as variables.
x1dc = x1d(2:end-1); 
if nVoffset==1
    [X,Y] = meshgrid(x1d, x1d);
elseif nVoffset==-1
    [X,Y] = meshgrid(x1dc, x1dc);
else
    error('Unsupported value of nVoffset!');
end

% some RHS function of the test in page 64 of A Multigrid Tutorial
sinX = sin(X);
sinY = sin(Y);
RHS = -3*sinX.*sinY;
rhs = RHS(:);

% the depth of v-cycle. Again, the number of cells has to be a power of 2.
q=log2(n);
% n is now the number of variables.
n=n+nVoffset;

% scale the RHS with h^2 and assign h=1 for the Laplacian operator.
rhs = h*h*rhs;
% assemble the LaplacianOperator
L2d = LaplacianOperator(n, 2, 1);
L2d = speye(size(L2d))-L2d;
% change the LaplacianOperator for the fact that
%  we know the values of the unknowns at the boundaries.
change2Identity = 0;
if change2Identity==1 && nVoffset==+1
    % the first column
    L2d(1:n,1:n) = speye(n);
    L2d(1:n,n+1:end) = speye(n,n*n-n);
    % the upper and lower boundary
    for i=2:n-1
        vid = (i-1)*n+1;
        L2d(vid,:) = sparse(1,n*n);
        L2d(vid,vid) = 1;
        vid = i*n;
        L2d(vid,:) = sparse(1,n*n);
        L2d(vid,vid) = 1;
    end
    % the last column
    rid = n*n-n+1:n*n;
    L2d(rid,:) = sparse(n,n*n);
    L2d(rid,rid) = speye(n);
end

% the exact solution we expect.
%sln = L2d\rhs;
% set the gobal matrices and initialize the MG 
mg_matrix = mg_matrix_init(L2d, q, 'p1', nVoffset);

% use zero as the initial guess.
uS = zeros(n*n,1); 
% the initial residual
res0 = norm(rhs,inf);
% the initial error
%err0 = norm(sln,inf);
for jj=1:maxIt
    if cycleType=='v'
%        uS = vcycle_block(uS, rhs, q);
        uS = vcycle(uS, rhs, q);
    elseif cycleType=='f'
%        uS = fcycle_block(uS, rhs, q);
        uS = fcycle(uS, rhs, q);
    else
        error('unknown cycle type!');
    end

    tmp = rhs-L2d*uS;
    res = norm(tmp,inf);
    resRel = res/res0;
 %   errRel = norm(uS-sln, inf)/err0;
    info = ['  iter #', int2str(jj), ': max-norm of Residual=', ...
        num2str(resRel)];
%     info = ['  iter #', int2str(jj), ': max-norm of Residual=', ...
%         num2str(resRel), '; max-norm of error=', num2str(errRel)];
    disp(info);
%    if (errRel<TOL)
    if (resRel<TOL)
        break
    end
end
nIters = jj;
errRel = 0;

end
